#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#include "MSTK.h"


#ifdef __cplusplus
extern "C" {
#endif

int MESH_Vol_Partition(Mesh_ptr mesh, Mesh_ptr *submeshes, int num) {
  int i, j;
  int nv, nf, nfv;
  MVertex_ptr mv, vertex;
  MFace_ptr mf;
  List_ptr mfverts;
  
  /* basic mesh information */
  nv = MESH_Num_Vertices(mesh);
  RepType type = MESH_RepType(mesh);
  nf = MESH_Num_Faces(mesh);

  /* Partition the mesh */
  int *part = (int *)MSTK_malloc(nf*sizeof(int));
  int ok = MESH_PartitionWithMetis(mesh, num, &part);
  /* collect entities of the submesh */
  int part_no;
  /* global to local vertex ID mapping */
  int *v_g2l = (int *)MSTK_malloc(nv*sizeof(int));
  /* for each partition */
  for(part_no = 0; part_no < num; part_no++) {
    /* clear to 0 */
    for( j = 0; j < nv; j++)
      v_g2l[j] = 0;
    int index_vertex = 0, index_face = 0;
    /* List of faces and vertices of this partition */
    List_ptr part_vertices=List_New(10), part_faces=List_New(10);
    int mkvid = MSTK_GetMarker();
    int max_nfv = 0;
    for (i = 0; i < nf; i++) {
      /* add face */
       if(part[i] == part_no) {
	 index_face++;
	 mf = MESH_Face(mesh,i);
	 List_Add(part_faces,mf);
	 mfverts = MF_Vertices(mf,1,0);
	 nfv = List_Num_Entries(mfverts);
	 if(max_nfv < nfv)
	   max_nfv = nfv;
	 /* add vertices */
	 for(j = 0; j < nfv; j++) {
	   mv = List_Entry(mfverts,j);
	   if (!MEnt_IsMarked(mv,mkvid)) {
	     v_g2l[MEnt_ID(mv)-1] = index_vertex++;
	     List_Add(part_vertices,mv);
	     MEnt_Mark(mv,mkvid);
	   }
	 }
	 List_Delete(mfverts);
       }
     }
    /* clear and release the marker */
    List_Unmark(part_vertices,mkvid); 
    MSTK_FreeMarker(mkvid);

    /* collect the coordinates */
    typedef double (*COORDS)[3];
    COORDS xyz =(COORDS) MSTK_malloc(index_vertex*3*sizeof(double));
    for (i = 0; i < index_vertex; i++) {
      vertex = List_Entry(part_vertices,i);
      MV_Coords(vertex,xyz[i]);
    }
    /* collect the face-vertex connectivity */ 
    int *nfvs = (int *)MSTK_malloc(nf*sizeof(int)); 
    int *fvids = (int *)MSTK_malloc(index_face*max_nfv*sizeof(int));
    for (i = 0; i < index_face; i++) {
      mf = List_Entry(part_faces,i);
      mfverts = MF_Vertices(mf,1,0);
      nfv = List_Num_Entries(mfverts);
      nfvs[i] = nfv;
      for (j = 0; j< nfv; j++) {
	mv = List_Entry(mfverts,j);
	fvids[i*max_nfv+j] = v_g2l[MEnt_ID(mv)-1];
	assert( fvids[i*max_nfv+j]<nv );
      }
    }
    List_Delete(part_vertices);
    List_Delete(part_faces);
    /* build the submesh */
    Mesh_ptr submesh = MESH_New(type);
    MESH_InitFromGenDesc(submesh,index_vertex,xyz,index_face,nfvs,&fvids,0,NULL,NULL,NULL,NULL);
    submeshes[part_no] = submesh;     
    MSTK_free(xyz);
    MSTK_free(nfvs);
    MSTK_free(fvids);
  }
  MSTK_free(v_g2l);
  MSTK_free(part);
  
  return ok;
}
  
#ifdef __cplusplus
}
#endif

